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Abstract 

Automated algorithms for derivation of amplitude equations in the vicinity of 
monotonia and Hopf bifurcation manifolds are presented. The implementation is 
based on Mathematica programming, and is illustrated by several examples. 



1 Introduction 



Amplitude equations describe behavior of evolutionary systems in the vicinity of thresh- 
old, or critical values of parameters that mark a point of qualitative change in behavior, 
like a transition from a quiescent state to convection, or from rigid mechanical equi- 
librium to vibration, or from a homogenious to patterned state of a material medium. 
Near the criticality, the behavior of systems of different physical origin is described by 
representative equations belonging to one of universal classes that are determined by the 
character of the transition. 

There is a number of celebrated "universal" equations which have been extensively 
studied in physical literature. Some of them were first suggested as models, others were 
derived as rational approximations in a certain physical context, and then discovered 
again in apparently unrelated problems. The general structure of these equations can be 
often predicted by symmetry considerations. The algorithms of their derivation, known as 
center manz/b/c? methods in mathematical literature [Guckenheimer & Holmes, 1983; Hale 
& Kogak, 1991], or elimination of "slaved" variables in physical literature [Haken, 1987], 
are well established, at least in standard cases. However, actual computations required 
for establishing quantitative relations between the coefficients of amplitude equations and 
parameters of underlying physical systems are repetitive and time-consuming, and have 
been carried out in scientific literature on a case-to-case basis. 

It is important to note that, except simplest cases when coefficients can be removed 
by rcscaling, different qualitative features of the behavior are observed in different param- 
eter ranges. When relations between coefficients of universal equations and measurable 
physical parameters are known, the domains determined by studying universal equations 
can be mapped upon the actual parameter space, and then further continued into a region 
where universal equations are no longer valid quantitatively but still faithfully describe 
qualitative features of behavior of a complex nonlinear system. 

The aim of this communication is to describe the Mathematica-hased automated al- 
gorithm for derivation of amplitude equations in either symbolic or numerical form. The 
algorithm makes use of Mathematica capabilities as a programming language [Wolfram, 
1991] that allow to define special functions carrying out complicated tasks in response to 
simple and well-defined inputs. The paper is organized as follows. The general algorithm 
based on multiscale bifurcation expansion is outlined in Section 2. The three functions 
designed to implement the algorithm of bifurcation expansion - Bif urcationCondition, 
BifurcationTheory, and Unfolding are described in Section 3. Several examples of us- 
age of these functions for derivation of amplitude equations for dynamical systems are 
given in Section 4. Further generalization of the algorithm to distributed systems will be 
described in forthcoming communications. 

2 General Algorithm 

2.1 Multiscale Expansion 

The general standard algorithm for derivation of amplitude equations is based on multi- 
scale expansion of an underlying dynamical system in the vicinity of a bifurcation point. 
Consider a set of first-order ordinary differential equations in R": 



du/dt — f (u, p). 



(1) 



where f (u; p) is a real- valued vector- function of an array of dynamic variables u, that is 
also dependent on an array of parameters p. We expand both variables and parameters 
in powers of a dummy small parameter e: 

u = uo + eui + e^U2 + . . . , 

p = po + epi + e^p2 + . . . . (2) 

We shall also introduce a hierarchy of time scales rescaled by the factor e^, thus replac- 
ing the function u(t) by the function of an array of rescaled time variables. Accordingly, 
the time derivative is expanded in a series of partial derivatives dk = d/dtk'. 

± = do + ed, + e'd2 + .... (3) 
at 

Let u = uo(po) be an equilibrium (a fixed point) of the dynamical system (|1]), i.e. a 
zero of the vector-function f (u; p) at a point p = po in the parametric space. The function 
f(u;p) can be expanded in the vicinity of uo,po in Taylor series in both variables and 
parameters: 

1 1 

f(u;p) = fu(u - Uo) + -fuu(u - Uo)^ + -fuuu(u - Uo)^ + . . . 

2 o 

+ fp(p - Po) + fup(u - uo)(p - Po) + -fuup(u - uo)^(p - Po) + . . . . (4) 

The derivatives with respect to both variables and parameters fu = df/du, = 
d'^f/du^, fp = df/dp, fup = d'^i/dudp, etc., are evaluated at u = uq, p = Pq. 
Using Eqs. (H) - (|) in Eq. ([l|) yields in the first order 

^oUi = fuUi + fpp^ 

The homogeneous linear equation 

Cu, = (f„ - d/dto)u, = (6) 

obtained by setting in Eq. © Pi = governs stability of the stationary state u = uo to 
infinitesimal perturbations. The state is stable if all eigenvalues of the Jacobi matrix fu 
have negative real parts. This can be checked with the help of the standard Mathematica 
function Eigenvalues. Computation of all eigenvalues is, however, superfluous, since 
stability of a fixed point is determined by the location in the complex plane of a leading 
eigenvalue, i.e. that with the largest real part. Generically, the real part of the leading 
eigenvalue vanishes on a codimension one subspace of the parametric space called a 
bifurcation manifold. The two types of codimension one bifurcations that can be located 
by local (linear) analysis are a monotonic bifurcation where a real leading eigenvalue 
vanishes, and a i/op/ bifurcation where a leading pair of complex conjugate eigenvalues 
is purely imaginary. Additional conditions may define bifurcation manifolds of higher 
codimension. 



(5) 



2.2 Monotonic Bifurcation 

A monotonic bifurcation is also a bifurcation of equilibria of the dynamical system (|l|). 
Setting ui = const, i.e. doUi = 0, allows to find a shift of the stationary solution due to 
small variations of parameters 

ui = -fu^fpp^. (7) 



The n X m array fu~ fp is recognized as a parametric sensitivity matrix; (m denotes the 
dimension of the parametric space). Continuous dependence on parameters can be used 
to construct a branch of equihbria which terminates at a bifurcation point po- 

On a monotonic bifurcation manifold, the matrix fu has no inverse. This means 
that one can neither construct a stationary solution at values of parameters close to this 
point, nor characterize the stability of the equilibrium in the linear approximation. The 
dynamics in the vicinity of the bifurcation manifold is governed by a nonlinear amplitude 
equation to be obtained in higher orders of the expansion. 

Generically, the zero eigenvalue is non-degenerate. Let U be the corresponding eigen- 
vector satisfying fuU = 0. Then 

ui =a(ti,t2,...)U (8) 

is the solution of Eq. that remains stationary on the fast time scale to- The amplitude 
a is so far indeterminate, and may depend on slower time scales. 

The inhomogeneous equation Eq. has solutions constant on the rapid time scale, 
provided its inhomogeneity does not project on the eigenvector U. This condition is 

K, = Utfpp^ = 0, (9) 

where is the eigenvector of the transposed matrix fu^ = Transpose [fu] satisfying 
fn^Ut = 0; we assume that the eigenvector is normalized: U^U = 1. Eq. @ defines the 
tangent hyperplane to the bifurcation manifold at the point P2 = Po- 

Before writing up the second-order equation, we require that the second-order devi- 
ation U2 remain constant on the rapid time scale (otherwise it may outgrow ui at long 
times). The dependence on slower time scales must be expressed exclusively through the 
time dependence of the amplitude a. Using Eq. (||), we write the second-order equation 
as ^ 

fuU2 = dialJ - fpP2 - afupUp^ - -a^fuuUU. (10) 
The solvability condition of this equation is 

dia = K2 + Aio — /ioci^ 

The parameters of Eq. (|1T]) are 

«:2 = U^fpP2, Ai = U^fupUp„ /io = -^U^UUU. (12) 

The indices correspond to the scaling of respective parametric deviations from the bifur- 
cation point. In a generic case, one can consider only parametric deviations transverse 
to the bifurcation manifold, and set pi = to satisfy Eq. (^; then Ai = 0. Thus, the 
generic equation for slow dynamics near the bifurcation manifold is 

dia = K2 — fioo^. (13) 

On the one side of the bifurcation manifold, where K2 has the same sign as fiQ, there are 
two stationary states a = ±^J K2/ /io- When viewed as a solution of Eq. ([13|), one of them is 
stable, and the other, unstable. The stable solution corresponds to a stable equilibrium 
of the original system Eq. (0), provided the rest of eigenvalues of the matrix fu have 
negative real parts. On the other side of the bifurcation manifold, where the signs K2 and 



(11) 



/io are opposite, there are no stationary states. The trajectory of the dynamical system 
is deflected then to another attractor, far removed from uq. Thus, the system undergoes 
a first-order phase transition when the bifurcation manifold is crossed. If the value of 
some dynamic variable or other characteristic of the solution is drawn as a function of 
parameters, the bifuration locus can be seen as the projection of a fold of the solution 
manifold on the parametric plane; accordingly, the generic monotonic bifurcation is also 
called a fold bifurcation. 

It may happen that the matrix fp vanishes identically. This would be the case when 
Uq is a "trivial" solution that remains constant at all values of parameters. Then Eq. (y) 
is satisfied identically, and Eq. ([TT|) reduces to 

dia = Aia — /ioa^- (14) 

This equation has two solutions, a = and a = Ai//i2, on both sides of the bifurcation 
manifold, but the two solutions interchange stability when this manifold is crossed. 



2.3 Unfolding of Higher-Order Bifurcations 

li Ho = 0, the expansion should be continued to the next order. The coefficient /Xq 
may vanish identically due to the symmetry of the original problem to inversion of u. 
Otherwise, it may be equal to zero at certain values of the parameters of the problem. 
Generally, the two conditions, Det[fu] = and /iq = are satisfied simultaneously on a 
codimension two manifold in the parametric space that corresponds to a cusp singularity. 

In order to continue the expansion, we restrict parametric deviations in such a way 
that the dependence on ti be suppressed. Deviations transverse to the bifurcation man- 
ifold have to be restricted by the condition K2 = 0, which is stronger than Eq. (|^). 
First-order parametric deviations pi parallel to the bifurcation manifold, which are still 
allowed by Eq. (^), should be now restricted by the condition Ai = 0. If the array 
p contains two parameters only, the conditions Ai = and Ki = imply, in a non- 
degenerate case, that first-order parametric deviations should vanish identically. When 
more parameters are available, parametric deviations satisfying both these conditions are 
superfiuous, since they correspond just to gliding into a closer vicinity of another point 
on the cusp bifurcation manifold in a higher-dimensional parametric space. Further on, 
we shall set therefore pi to zero identically. 

The dynamics unfolding on a still slower time scale ^2 should be determined from the 
third-order equation 

fuUg = 52aU - fpPg - afupUp2 - afuuUu2 - ^a^fuuuUUU. (15) 

The second-order function U2 has to be found by solving Eq. (p!OD, now reduced to the 
form 

fuU2 = fpP2 + ^a^f^^UU. (16) 

Only the solution of the inhomogeneous equation, which does not project on the eigen- 
vector U, is relevant. It can be expressed as 

U2 = + a2u(°\ (17) 

The solvability condition of Eq. ([T5|) is obtained then in the form 

= K3 + A2a - z/qo^, (18) 



where 



K3 = U^fpPg, 

A2 = Utf,pUp2 + Utf,,UU^'\ 

z/0 = -^Utf,,,UUU-UtUUUf . (19) 
o 

Eq. ([TSD presents the parametric unfolding of dynamics in the vicinity of a cusp 
bifurcation. Three equihbria - two stable and one unstable - exist in the cusped region 

^2>0, \^3\<^^j^- (20) 

Outside this region, there is a unique stable equilibrium. A second- order phase transition 
occurs when the parameters change in such a way that A2 crosses zero. Other transitions 
occuring in the vicinity of the cusp bifurcation are weakly first-order. 

The condition uq = defines a singular bifurcation manifold of codimension three. 
Again, it is possible to fix parametric deviations to suppress the dynamics on the scale 
t2, and obtain in the next order a quatric equation that represents the unfolding of the 
butterfly singularity. The procedure can be continued further if a sufficient number of 
free parameters is available. The unfolding of a codimension q singularity is presented 
by a polynomial of order q + 1: 

q-l 

dqtt = (^q^p+ia^ - o-oa''"*"\ (21) 

p=0 

where the parameters cxfc depend on parametric deviations proportional to e*^. 
2.4 Hopf Bifurcation 

At the Hopf bifurcation point, the parametric dependence of the stationary solution 
u = Uo(p) remains smooth; linear corrections can be obtained from the stationary eq. (|^), 
and higher corrections from higher orders of the regular expansion (|). In order to simplify 
derivations, we shall eliminate this trivial parametric dependence by transforming to a 
new variable u = u — uo(p). The resulting dynamic system, dix/dt = f(u, p) has the 
same form as (||) but contains a modified vector-function f(u, p) = f(u + uo,p). Since, 
by definition, uq satisfies f(uo,p) = 0, u = is a zero of f(u, p). Now we can drop the 
hats over the symbols and revert to the original form (|I|) while keeping in mind that 
u = is a stationary solution for all p and, consequently, all derivatives fp, fpp, etc., 
computed at u = vanish. 

At a Hopf bifurcation point the Jacobi matrix fu has a pair of imaginary eigenvalues 
A = iituo- The first order equation (||) has a nontrivial oscillatory solution 

ui = a(ti, t2, . . -Mto) + CO.; $(to) = e^-°*«U (22) 

with an arbitrary complex amplitude a{ti,t2, ■ ■ .) dependent on slower time scales tj, i > 
0; U is the eigenvector of fu with eigenvalue iuo- 



fuU = ^CUqU. 



(23) 



The function $(to) and its complex conjugate $*(to) a^re two eigenfunctions of the hnear 
operator £ with zero eigenvalue. The operator £ is defined in the space of 27r/co'o-periodic 
complex-valued vector-functions with the scalar product defined as 

(u,v) = -^/ u\t).w{t)dt. (24) 
Ztx Jo 

The eigenfunctions of the adjoint operator = i^^ + d/dto are 

$t(to) = e-^'^«*oU^ 

and its complex conjugate; is the eigenvector of fu^ with the eigenvalue iuo. 
The second-order equation can be written in the form 

£u2 = diui - ^fuuUiUi - fupUipi (25) 

The inhomogeinity of this equation contains both the principal harmonic, e^'^°^°, con- 
tributed by the linear terms diUi and fupUiPi, and terms with zero and double frequency 
coming from the quadratic term ifuuUiUi. The scalar products of the latter terms with 
the eigenfunctions $(to), '^'*(^o) vanish, and the solvability condition of eq. ( pS]) is ob- 
tained in the form: 

dia = Aia, (26) 

where Ai is given in Eq. (|12D (note that this parameter is now complex). 

Equation (|26|) has a nontrivial stationary solution only if the real part of Ai vanishes. 
This condition defines a hyperplane in the parametric space tangential to the Hopf bifur- 
cation manifold. The imaginary part of Ai gives a frequency shift along the bifurcation 
manifold. In order to eliminate the tangential shift, we set, as before, pi = 0. Then 
eq. (p6D reduces to dia = 0, so that the amplitude can evolve only on a still slower scale 
^2- The second-order function U2 has to be found by solving Eq. (P^D, now reduced to 
the form 

£u2 = -^(a'e''"o*«fuuUU+ | a p f^uUU* + c.c). (27) 
The solution of this equation is 

U2 = « r fu^'fuuUU* + a'e^'^^'^ii^ - 2^a;o)"'fuuUU + c.c.]. (28) 
In the third order, 

1 

-1 
6 

The amplitude equation is obtained as the solvability condition of this equation. Only the 
part of the inhomogeinity containing the principal harmonic contributes to the solvability 
condition, which takes the form 

820, = — z/q I a 1^ a, (30) 

where A2, t'o are defined as 

A2 = Utf„pUp2, 

^0 = -^Utf,,,UUU* + Utf,,U(f,-^f,,UU*) + 

iutf,,U*((f„ - 2tujo)-%^UU). (31) 

The expansion has to be continued to higher orders, after readjusting the scaling of 
parametric deviations, if the real part of u vanishes. 



£U3 = d2Ui - fupUiP2 - fuuUlU2 - -fuuuUlUiUi. (29) 



3 Automated Generation of Amplitude Equations 



3.1 Localization of Bifurcation Points 

Derivation of an amplitude equation for a dynamical system should be preceded by local- 
ization of bifurcation manifolds. The function Bif urcationCondition [f , u, options] 
generates a set of algebraic equations defining a bifurcation point of a specified type 
(monotonic or Hopf). Bif urcationCondition has two arguments: f denotes an array 
of functions in Eq. (|1|), and u denotes an array of variables. 

The type of a bifurcation point is specified by the option Special with the default 
value Special -> None that corresponds to a codimension one monotonic bifurcation. 
The set of equations defining the bifurcation manifold is produced then by adding the 
condition of vanishing Jacobian Det{di/du) = to the stationarity conditions f(u) = 0. 

If the dynamical system is one-dimensional, the option Special can be set to an inte- 
ger q; then the set of equations defining monotonic bifurcations of a higher codimension is 
produced by equating to zero subsequent derivatives of the function f{u) up to the order 
q: (/ = 0, df/du = 0, d^f/dv? = 0, . . . , d'^f/du'^ = 0). For n > 1 this option is not im- 
plemented, since the required condition depend on the eigenvectors of the Jacobi matrix 
at the bifurcation point. The symbolic form of the conditions of higher codimension bi- 
furcations can be obtained in this case with the help of the function Bif urcationTheory 
(see below). 

Specifying Special -> Hopf produces the condition of the Hopf bifurcation. For an 
n-dimensional dynamical system, the additional condition supplementing the stationar- 
ity conditions is vanishing of the {n — l)-th Hurwitz determinant. This is a necessary 
condition of Hopf bifurcation, which is also sufficient when there are no eigenvalues with 
positive real parts; in the latter case, the Jacobian must be positive. This case is usually 
the only interesting one, as it implies that the stationary state in question is indeed stable 
at one side of the Hopf bifurcation manifold. A more complex algorithm [Guckenheimer, 
Myers & Sturmfels, 1996] can be implemented to obtain a more general necessary and 
sufficient condition for existence of a pair of imaginary eigenvalues. 

3.2 Function Bif urcationTheory 

The procedure of derivation of a sequence of amplitude equations described in the pre- 
ceding Section is implemented by the function Bif urcationTheory. 

Bif urcationTheory [f , u ,p, t, amp, {U,U^}, order, options] constructs in a 
symbolic form a set of amplitude equations on various time scales for a nonlinear dy- 
namical system du/dt = f(u, p) in the vicinity of a bifurcation point, according to the 
algorithm described in Section 2. 

Bif urcationTheory has the following arguments: 

• f denotes an array of functions in eq. (|l|); 

• u denotes an array of variables; 

• p denotes an array of bifurcation parameters; 

• t denotes the time variable; 

• amp is a name of the amplitude function to appear in the amplitude equations; 



• U and are names for the eigenvectors of a hnearized problem and its adjoint, 
respectively; 

• order is the highest order of the expansion; if it is higher than the codimension 
of the bifurcation manifold, corrections to the principal (lowest order) amplitude 
equation are produced. 

The function BifurcationTheory admits two options. The option Special specifies 
the type of the bifurcation: Special -> None (default) corresponds to a monotonic bi- 
furcation, and Special -> Hopf, to a Hopf bifurcation. The option BifurcationLocus, 
if set to True, causes an intermediate printout of equations determining the location of 
a bifurcation point of codimension order. The default value of this option is False. 

BifurcationTheory needs not to be implemented in actual computations, but only 
serves to display the internal algorithm and to generate symbolic conditions for higher 
codimension bifurcations. The output of this function, that can be very bulky, contains 
all applicable formulae from Section 2 written in Mathematica format. 

In the following example, BifurcationTheory is applied to a codimension two (cusp) 
monotonic bifurcation: 

In[lj~ 

BifurcationTheory [f , u , p , t , amp , {U , Ut } , 2 , 
BifurcationLocus -> True] 

{f[u, p] ==0, Det[f(0'i) [u, p]] ==0, 
Ut . f(2'0)[u, p] . U . U 

== 0} 

2 

Outll]^ 

{0 == f [u, p] , 

{0 == ut . f(0'i)[u, p] . p[l] / Ut . U, {}, t[0]}, 
{amp(i'0)[t[l] , t[2]] == Ut . f(0'i)[u, p] . p[2] / Ut . U + 
Ut . f(°'2)[u, p] . pCl] . p[l] / (2 Ut . U) + 
(ampEtEl], t[2]] Ut. f(^'i)[u, p] . U . p[l]) / Ut . U + 
(amp[t[l], t[2]]2 Ut . f(2.o)[u, p] . U . U) / (2 Ut . U) , 
{u[2][t[l], t[2]] -> 

InverseOperator[f(i'°) [u, p]] . (-f^^'^Eu, p] . p[2] - 

f(0'2)[u, p] . pCl] . p[l] / 2 - 

amp[t[l], t[2]] f^i)[u, p] . U . p[l] - 

(ainp[t[l], t[2]]2 f(2.o)[u, p] . U . U) / 2 + 

U amp(i'O) [t [1] , t[2]])},t[l]}, 

{amp(O'i) [t [1] , t[2]] == 

-(Ut . u[2](i'0)[t[l] , t[2]] / Ut . U) + 

Ut . f(°'^)[u, p] . p[3] / Ut . U + 

Ut . f(°'2)[u, p] . p[l] . p[2] / Ut . U + 

(ampEtEl], t[2]] Ut . f(i'i)[u, p] . U . p[2]) / Ut . U + 

Ut . f(i'i)[u, p] . u[2][t[l], t[2]] . p[l] / Ut . U + 

(amp[t[l], t[2]] Ut .f(2.o)[u, p] .U . u[2] [t [1] ,t [2] ] ) / Ut.U + 

Ut . f^-'^'O^Cu, p] . U . U . U / (6 Ut . U) + 

(amp[t[l], t[2]] Ut .f(i'2)[u, p] .U .p[l] .p[l]) / (2 Ut . U) + 



(amp[t[l], t[2]]2 Ut .f(2'i)[u, p] .U . U . p [1] ) / (2 Ut . U) + 
(amp[t[l], t[2]]3 Ut .f(3.o)[u, p] .U . U . U) / (6 Ut . U)}, 
f(i'0)[u, p] . U == 0, 
Transpose [f(i'°) [u, p]] . Ut == 0}} 

The Mathematica output shown here starts with an intermediate printout (generated 
by setting the option Bif urcationLocus -> True) that contains the set of equations 
determining the bifurcation point. The output proper is a nested array containing the 
equations in the consecutive orders of the expansion, their solvabihty conditions and 
solutions. 

The first element (in the zero order) defines the fixed point of the system. The next 
element (in the first order) is an array that starts with the condition (P) defining the 
tangent hyperplane to the bifurcation manifold. The solution in this order, proportional 
to the eigenvector U is not presented, and replaced by an empty array. The last element of 
the first-order array is the time scale t [0] . The structure of the consecutive (up to order 
next to last) arrays is similar. The first element in each order is the amplitude equation 
on the corresponding time scale; the second element is the solution for the deviation 
u[i] from the fixed point, and the last is the time scale t[i-l]. The last array of the 
output, corresponding to the highest order of the expansion, contains only the amplitude 
equation and the equations determining the eigenvectors U,Ut. 

The same function generates a similar (but longer) output for case of Hopf bifurcation 
if the option value Special -> Hopf is specified. The output is shortened by specifying 
automatically some parametric conditions, say p[l] -> 0. 

3.3 Function Unfolding 

The actual derivation of the amplitude equation is carried out by the function 
Unfolding[f, u,uO, p, pO, amp, order, options]. The arguments of Unfolding 
are: 

• f denotes an array of functions in eq. (|I|); 

• u denotes an array of variables; 

• uO = uq denotes an array of their values at the bifurcation point; 

• p denotes an array of bifurcation parameters; 

• pO = po denotes an array of their values at the bifurcation point; 

• amp is a name of the amplitude function to appear in the amplitude equations; 

• order = q denotes the highest order of the expansion. 

The only option of this function. Special, has the same values as for 
Bif urcationTheory. 

The additional arguments of Unfolding, as compared to Bif urcationTheory, are the 
bifurcation values of the variables and parameters that can be given either sympolically 
or numerically. The function first tests whether the provided bifurcation point is indeed 
an equlibrium of the dynamical system. Then it computes the Jacobi matrix 3 = di/ du 
and checks whether the conditions of either monotonic or Hopf bifurcation are verified. 



If the test is passed, the eigenvectors U, of the Jacobi matrix and its transpose are 
computed, and the rest of computations leading to an amphtude equation in the desired 
order are carried out using Bif urcationTheory as a subroutine. 



4 Examples 



In this section we present several examples of automated derivation of amplitude equa- 
tions for particular dynamical systems (some additional examples are found in [Pismen, 
Rubinstein & Velarde, 1996]). 

4.1 Monotonic Bifurcations in the Lorenz Model 

The well-known Lorenz system [Lorenz, 1963], that has been originally suggested as a 
qualitative model of cellular convection, exhibits rich dynamic behavior including periodic 
and chaotic motion. It also possesses a monotonic bifurcation that corresponds to a 
primary transition from the quiescent state to convection. In usual notation, the Lorenz 
system is written as 



We start with defining of the array f containing the right-hand sides of the above equa- 
tions: 



Lorenz = {-sigma x + sigma y, Rx-xz-y, xy-bz}; 

and call the function Bif urcationCondition. As the option Special is not specified, a 
monotonic bifurcation is computed by default: 



We choose R as the bifurcation parameter, and solve the above equations with respect 
to X, y, z, R: 

In[4]:^ 

bpts = Solve [be, { x, y, z, R}] 

Outf4j= 

{{R -> 1, z -> 0, y -> 0, x -> 0}, 
(R -> 1, z -> 0, y -> 0, X -> 0}, 
{R -> 1, z -> 0, y -> 0, X -> 0}} 

The solution is degenerate, so there is in fact a single point of monotonic bifurcation. 
Now the function Unfolding is called to derive the amplitude equation. 



dx/dt 
dy/dt 
dz/dt 



-a{x - y), 
x{R - z)-y, 
—bz + xy. 



In[2]: 




In[5]:^ 

bp = First [bpts] ; 

Unfolding [Lorenz,{x,y,z},{x,y,z} /. bp, 
{R},{R} /. bp,amp[t] ,2] 

Unfolding: :parm: 

The list of parameters {R} 

has a unsufficient length for the bifurcation point of 
codimension 2. 

Out[5]^ 

{{x == amp [t [2] ] , y == amp [t [2] ] , z == 0} , 
{amp'[t[2]] == - Sigma amp[t [2]]3/(b (1 + 1/sigma)) + 
amp[t[2]] R[2] /(b (1 + 1/sigma))}} 

The output is preceded by a warning message, because the number of specified bifur- 
cation parameters is smaller the codimension of the specified bifurcation point. The 
amplitude equation is obtained nevertheless, but the user is warned that the unfolding is 
not complete. In a generic case, this might have caused a major trouble, since it would 
be impossible to satisfy the parametric condition K2 — 0. The Lorenz system possesses, 
however, an inversion symmetry that causes this condition to be verified identically, and 
therefore a single bifurcation parameter is in fact sufficient. 

The output is a nested array. The first clement is a list expressing the original 
dynamical variables through the amplitude (to the first order in the expansion parameter 
e). The second element contains the desired amplitude equation. The actual scafing of 
the amplitude in terms of parametric deviations from the bifurcation point is indicated by 
the presence in the amplitude equation of the second-order deviation R [2] . The result 
proves that the bifurcation is always supercritical at physical (positive) values of the 
parameters. 

It is possible to find a higher-order correction to the amplitude equation derived above 
by specifying the maximal order of the expansion equal to four: 

Inl6]:= 

Unfolding [Lorenz, {x,y,z},{x,y,z} /. bp, 
{R},{R} /. bp,amp[t] ,4] 

Unfolding: :parm: 

The list of parameters {R} 

has a unsufficient length for the bifurcation point of 
codimension 2. 

Out[6]= 

{{x == amp[t[2], t[4]], y == amp[t[2], t[4]], z == 0}, 
{amp(i'O) [t [2] , t[4]] == 

- Sigma amp[t[2], t[4]]^ / (b (1 + 1 / sigma)) + 
amp[t[2], t[4]] R[2] / (b (1 + 1 / sigma)), 

amp(O'i) [t [2] , t[4]] == 

(Sigma (b - 2 Sigma -2b sigma - 2 sigma^) 
amp[t[2], t[4]]5) / (b^ (1 + sigma)^) + 
(sigma (-b + 2 sigma + 3 b sigma + 2 sigma^) 
amp[t[2], t[4]]3 R[2] / (b^ (1 + sigma)^) - 



(sigma^ amp[t[2], t[4]] R[2]2) / (1 + sigma)^}} 



4.2 Hopf Bifurcation in the Brusselator Model 

The next example is a well-known Brusselator model [Nicolis & Prigogine, 1977] of chem- 
ical oscillations: 

dz/dt = a — {1 + b)z + uz"^, 
du/dt = bz — uz^. 

The variables u, z denote concentrations of the "activator" and "inhibitor" species. We 
define the array brusselator containing the r.h.s. of the system: 

In[7j:^ 

brusselator = {a - (1 + b) z + u z"2,b z - u z^2}; 

This system always has a unique stationary solution z — a, u — b/a. The condition for 
Hopf bifurcation point is defined as follows: 

In[8]:= 

bpt = First [ Solve [ Bif urcationCondition[ brusselator, {z,u}. 
Special -> Hopf], {z,u,b}]] 

Out[8]= 

(b -> 1 + a^, u -> (1 + a^) / a, z -> a} 

Now the function Unfolding with the option Special -> Hopf is called to produce an 
amplitude equation: 

In[9j:^ 

Unfolding [brusselator, {z,u},{z,u} /. bpt, 
{b}, {b} /. bpt, ampLt], 2, Special -> Hopf] 

Out[9]^ 

{{z == a - a E^^*[°l ainp[t[2]] / (-1 + a) - 
a E-^^*[°l Conjugate [amp [t [2]]] / (I + a) , 
u == (1 + a2) / a + E^^*[°l amp[t[2]] + 
£-1 a t[o] Conjugate [amp [t [2] ] ] } , 

If[a2 > 0,{amp^ [t[2]] == (amp[t[2]] b[2]) / 2 + 
((4 - 6 I a - 7 a^ - 3 I a^ + 4 a^) amp[t[2]]2 
Conjugate [amp [t [2] ]] ) / 
(6 (-1 + I a) a (-1 + a))}]} 

The If conditional clause which seems to be redundant is generated because Mathematica 
understands all symbolic variables as complex quantities. The coefficient at the nonlinear 

term is generated in a somewhat clumsy form but it is easy to see that its real part 
simplifies to — {1 + ^a^) / {1 + a^) , and is negative definite, hence, the bifurcation is always 
supercritical. 

It is possible to find a fourth-order correction to the above equation: 
In[10]:= 

Unfolding [Brussel,{z,u},{z,u} /. bpt, 

{b}, {b} /. bpt, amp[t], 4, Special -> Hopf] 



Out[10]= 

{{z == a - a E^^^M amp[t[2], t[4]] / (-1 + a) - 
3^ E-iat[o] Conjugate [amp [t [2] , t[4]]] / (I + a), 
u == (1 + a2) / a + E^^^M amp[t[2], t[4]] + 
£-1 a t[o] Conjugate [amp [t [2] , t [4] ] ] } , 

lf\_3? > 0,{amp(^'°) [t[2] , t[4]] == (amp[t[2], t[4]] b[2])/2 + 

((4 - 6 I a - 7 a^ - 3 I a^ + 4 a'') amp[t[2], t[4]]2 

Conjugate [amp [t [2] , t[4]]]) / (6 (-1 + I a) a (-1 + a)), 

amp(O'i) [t [2] , t[4]] == - (I amp[t[2], t[4]] b[2]2) / (4 a) + 

((-28 + 54 I a + 25 a^ + 27 I a^ - 28 a^) 

amp[t[2], t[4]]2 b[2] Conjugate [amp [t [2] , t[4]]]) / 

(72 a2 (-1 + a) (I + a)) + 

((112 I + 576 a - 908 I a^ + 1155 I a"^ - 2088 a^ - 
1241 I a^ + 576 a^ + 112 I a^) amp[t[2], t[4]]3 
Conjugate [amp [t [2] , t[4]]]2) / 
(432 a^ (-1 + a)2 (I + a)^)}]} 

4.3 Exothermic Reaction in a Stirred Tank Reactor 

The following example involves somewhat heavier computations, and necessitates the use 
of implicit functions. Consider a dynamical system describing an exotermic reaction in 
a continous stirred tank reactor [Uppal, Ray & Poore, 1974; Pismen, 1986]: 

dx/dt = {l — x)e^ — mx, 
gdy/dt = h{l~x)ey -my. (32) 

The variables x, y denote the conversion and dimensionless temperature, respectively; h 
is the exotermicity parameter, m is the dimensionless flow rate, and g is the thermal 
capacitance factor that equals to unity in an adiabatic reactor and decreases with inten- 
sified cooling. This model exhibits both monotonic and Hopf bifurcations. Equilibria of 
Eq. (p2D do not depend on the parameter g, so conditions of a monotonic bifurcation 
may depend on two independent parameters m, h only; the maximal codimension is two 
(a cusp point). 

Monotonic bifurcation Consider first a bifurcation manifold of codimension 
one (a fold line). In order to localize this bifurcation we call the function 
Bif urcationCondition. 

In[ll]:= 

react = {(1-x) Exp[y] - m x,(h (1-x) Exp[y] -my) / g}; 
eqs = First [Bif urcationCondition [react /. {g ~> l},{x,y}]] 

Out[llJ= 

{E^ (1-x) - m X == 0, E^ h (1 - x) - m y == 0, 
E^ m - E^ h m + m^ + E^ h m X == 0} 

The bifurcation locus parametrized by the stationary value yg of the variable y can be 
found with the help of the standard Mathematica function Solve. This parametrization 
is advantageous, since if one of the actual parameters h or m was used, either two or 
no solutions were obtained in different parametric domains, and the form of the solution 
was more complicated. 



In[12]:= 

rl = Last [Solve [eqs , {x,ni,h}] ] // Simplify /. y -> ys 
Solve : : svars : 

Warning: Equations may not give solutions for all 

"solve" variables. 

Outfl2j= 

{m -> E^^ / (-1 + ys), h -> ys^ / (-1 + ys) , x -> 1 - 1 / ys} 

Solve generates a warning message but nevertheless produces a correct result. The 
admisible values of ys are restricted by the condition ys > I. 

The function Unfolding is now called to produce the amplitude equations on two 
time scales: 

In[13]:= 

rll = Join[rl,{y -> ys}] ; 

Unfolding [react /. {g -> 1} , {x,y} , {x,y} /. rll, 
{h,m},{h,m} /. rll, amp[t], 2] 

Out [13]= 

{{x == 1 - 1 / ys + (-1 + ys) ainp[t[l], t[2]] / ys^, 

y == ys + amp[t[l] , t[2]]}, 

{amp(i'°)[t[l] , t[2]] == 

E^^ (2 - ys) amp[t[l] , t[2]]2 / 

(2 (-1 + ys)) + E^^ h[2] - ys m[2] , 

amp(°'i) [t [1] , t[2]] == 

E^" (3-2 ys) ainp[t[l], t[2]]3 / (6 (-1 + ys)) + 
(amp[t[l], t[2]] m[2] (2 E^^ h[2] - 2 E^^ ys h[2] + 
E^^ ys2 h[2] - ys m[2])) / ys}} 

The second element of the output contains two amplitude equations on two time scales 
t[l] and t [2] . The first equation has the form (0), while the second equation gives 
a correction in the next order. Since this is a generic codimension one bifurcation, no 
relations between the parametric deviations m [2] and h [2] are generated, so that both 
can be varied independently. 

Cusp point The loci of fold bifurcation in the parametric plane h, m, drawn as a 
parametric plot, are shown in Fig. ||. The two branches join with a common tangent 
at the cusp point. In the vicinity of this point, the applicable amplitude equation will 
include a restriction on deviations of both parameters that specify the direction of this 
common tangent. In this simple example, the cusp point can be located analytically with 
the help of the function Bif urcationCondition using the fact that in the steady state 
the system (^) can be reduced to a singe equation by setting x = y/h: 

In[14]:= 

reactl = react /. {g -> l,x -> y/h} // Simplify 

Out[14]= 

{(E^ h - E^ y - m y) / h, E^ h - E^ y - m y} 

The set of equations defining the location of the cusp point is produced by setting the 
value of the option Special equal to the codimension of the bifurcation point: 



In[15]:= 

eqsc = BifurcationCondition [Last [reactl] , {y} , Special -> 2] 

Out[15]= 

{E^ h - y - m y == 0, - + h - m - y == 0, 
-2 E^ + E^ h - Ey y == 0} 

This system of equations is solved by combining the Mathematica functions Solve and 
Eliminate, yielding the values of the variable y and the two parameters at the cusp point 
(some irrelevant warning messages appear because of the common factor i?^): 

In[16]:= 

yc = First [Solve [Eliminate [eqsc ,{m,h}] ,y] ] ; 
bpc = Join[sl , First [Solve [eqsc /. yc,{m,h}]]] 

Solve : : ifun: 

Warning: Inverse functions are being used by Solve, so 
some solutions may not be found. 

Infinity: : indet : 

Indeterminate expression (-Infinity) encountered. 

Out[16j= 

{y -> 2, m -> E^, h -> 4} 

The corresponding value of the variable x at the cusp point is added, and the function 
Unfolding is called. 

In[17]:= 

bp = Join[{x -> y/h} /. bpc, bpc]; 

Unfolding [react /. {g -> l},{x,y},{x,y} /. bp, 

{h,m},{h,m} /. bp, ajnp[t], 2] 

Out[17j= 

{{x == 0.5 + 0.25 amp[t [2]] , y == 2 + 1 . ainp[t[2]]}, 
{amp'[t[2]] == -1.23151 amp[t[2]]3 + 7.9304 h[3] + 
1. amp[t[2]] m[2], h[2] -> 0.270671 m[2]}} 

Note that the output in this case contains not only the amplitude equation but also a 
relation between the second-order parametric deviations which is satisfied at the tangent 
line defined by Eq. (|). 

It is possible to find a correction to the above equation. To this end, the maximal 
order of expansion must be reset to three: 

In[18]:= 

Unfolding [react /. {g -> l},{x,y},{x,y} /. bp, 
{h,m},{h,m} /. bp, ajnp[t], 3] 

Out[18j= 

{{x == 0.5 + 0.25 amp[t[2], t[3]], 
y == 2 + 1 . amp [t [2] , t [3] ] } , 
{amp(^'°) [t [2] , t[3]] == 

-1.23151 ainp[t[2], t[3]]^ + 7.9304 h[3] + 
1 . amp [t [2] , t [3] ] m [2] , 



amp(°'i) [t [2] , t[3]] == 

-0.615755 ainp[t[2], t[3]]^ + 

7.65973 amp[t[2], t[3]] h[3] + 

1.5 ainp[t[2], t[3]]2 m[2] - 0.135335 m[2]2, 

h[2] -> 0.270671 m[2]}} 

Hopf bifurcation The dynamical system ( p2D can also undergo a Hopf bifurcation. 
First, we compute the Hopf bifurcation manifold. Taking into account the stationary 
relation x = y/h and applying the function Bif urcationCondition yields the set of 
equations defining the locus of Hopf bifurcation: 

In[19]:= 

eqs = Simplify [Rest [ 

Bif urcationCondition [react ,{x,y} , 

Special -> Hopf]] /. {x -> y / h}] /. {y -> yO} Out[19]= 

{(Ey° h - yO - m yO) / g == 0, 

(Ey° g + yO - h + m + m g) / g == 0} 

These equations are solved to express the values of the parameters m and h through the 
stationary value y = yo and the remaining parameter g = go; yo, go are thus chosen to 
parametrize the 2d Hopf bifurcation manifold in the 3d parametric space. 

In[20]:= 

bph = Simplify [First [Solve [eqs , {in,h}] ] ] 

Out[20j= 

{h -> yO + (g yO) / (-1 - g + yO) , 
m -> (Ey° g) / (-1 - g + yO)} 

The loci of Hopf bifurcation in the parametric plane h, m at several chosen values of g are 
shown in Fig. H Outside the cusped region, the unique stationary state suffers oscillatory 
instability within the loop of the Hopf curve (this is possible at g < ^ only). Within the 
cusped region, the solutions on the upper fold are unstable below the Hopf bifurcation 
line. The solutions at the lower fold may also suffer oscillatory instability at lower values 
of 5f. 

We write a set of rules defining the Hopf bifurcation manifold, and simplify the 
dynamical system in its vicinity: 

In[21]:= 

rO = {x -> yO/h,y -> yO,g -> gO} /. (bph /. g -> gO) ; 
fh = Simplify [(react / E"yO) /. bph] 

Out[21J= 

{(Ey - Ey X + (Ey° g x)/(i + g - yo)) / Ey°, 

(Ey° g y + E^ yO - E^ X yO - Ey yO^ + E^ x yO^) / 
(Ey° (g + g^ - g yO))} 

After this preparation, the function Unfolding is called to produce the amplitude equa- 
tion. 

In[22]:= 

Unfolding [fh,{x,y},{x,y} /. rO, 

{g}>{g} /■ rO, a[t], 2, Special -> Hopf] 



Out[22]= 

{{x == (-1 - gO + yO)/(-l + yO) - 

(£(I Sqrt[-1 + yO - gO yO] t[0])/(l + gO - yO) 

gO (1 - yO - I Sqrt[-1 + yO - gO yO] ) a[t[2]])/ (-yO + yO^) - 

Sqrt[-1 + yO - gO yO] t[0])/(l + gO - yO) 

gO (1 - yO + I Sqrt[-1 + yO - gO yO] ) 
Conjugate [a [t [2] ]])/(-yO + yO^) , 
y == yO + 

£(lSqrt[-l + yO-gOyO]t[0])/(l + gO-yO) a [t [2] ] + 

sqrt[-i + yo - go yo] t[o])/(i + go - yo) Conjugate [a [t [2] ] ] } , 
If [(-1 + yO - gO yO)/(l + gO - yO)^ > 0, 
{a'[t[2]] == 

((12 I - 24 I gO - 24 I yO + 40 I gO yO - 
4 I gO^ yO + 14 I yO^ - 21 I gO yO^ + 

4 I gO^ yO^ - 2 I yO^ + 2 I gO yO^ + 
12 Sqrt[-1 + yO - gO yO] - 

24 gO Sqrt[-1 + yO - gO yO] - 
15 yO Sqrt[-1 + yO - gO yO] + 

25 gO yO Sqrt[-1 + yO - gO yO] + 
2 gO^ yO Sqrt[-1 + yO - gO yO] + 

5 yO^ Sqrt[-1 + yO - gO yO] - 

8 gO yO^ Sqrt[-1 + yO - gO yO] ) a[t[2]]2 
Conjugate[a[t[2]]]) / (12 (1 - yO + gO yO) 
(-1 + I yO - I gO yO - Sqrt[-1 + yO - gO yO] + 
yO Sqrt[-1 + yO - gO yO])) + 
((2 - 3 yO + gO yO + yO^) 

(1 - yO - I Sqrt[-1 + yO - gO yO] ) a[t[2]] g[2] 
)/(2 (1 + gO - y0)2 

(1 - yO + gO yO - I Sqrt[-1 + yO - gO yO] + 
I yO Sqrt[-1 + yO - gO yO]))}]} 

This amplitude equation can be used to determine stability of a limit cycle with the 
amplitude a in the vicinity of the Hopf bifurcation. First, we take note that only a part 
of the parametric plane yo, go above the thick solid curve in Fig. || is actually available, 
since the condition of positive oscillation frequency requires 

yo - govo - 1 > 0. (33) 

Within this region, the limit cycle is stable (the bifurcation is supercritical) if the real 
part of the coefficient at the nonlinear term is negative. Extracting this coefficient from 
the above output and separating the real and imaginary parts brings this condition to 
the form 

3^0 + 2g^ + 2yo - igpyo - 2g^yo -yl + 2goy^ - 1 ^ ^ 

4(?/o - 1 - gQ){yo - goVo - i) 

Noting that the denominator of the above fraction is positive whenever the inequality ( P^D 
is verified, one can determine the stability boundary by equating the numerator in ( p4D 
to zero, and combining the result with (|33|) . The stability boundary in the plane {go, yo), 
separating the regions of subcritical and supercritical bifurcation is shown in Fig. Q the 
unphysical part of the curve dipping below the boundary of positive frequency is shown 



by a dashed line. The region of subcritical bifurcations consists of two disconnected 
parts, of which the lower one lies within the region of unique stationary states and on the 
lower fold of the solution manifold, and the upper one, on the upper fold in the region 
of multiple solutions. 

5 Conclusion 

Automatic algorithms provide basic tools for comprehensive nonlinear analysis in cases 
when the required procedures are known in principle but cannot be practically imple- 
mented because of cumbersome and heavy computations involved. The examples of 
analytic derivation of amplitude equations given above are relatively simple, but even 
they involve rather lengthy computations. The algorithm works much faster if all nec- 
essary data are given numerically; this would be the only possibility in more complex 
cases when analytical expressions for bifurcation loci are unavailable. Since, however, the 
underlying algorithm is symbolic, parametric deviations would appear in the resulting 
amplitude equations in a symbolic form, even when the rest of coeficients are numerical. 

In the subsequent communications we shall describe how the above algorithm can be 
generalized to be also applicable to distributed systems, and to carry out the following 
operations: 

• derivation of long-scale equations through dimensional reduction; 

• construction of loci of symmetry-breaking bifurcations; 

• construction of amplitude equations at stationary and oscillatory symmetry break- 
ing bifurcations; 

• detection and analysis of degenerate bifurcations. 

For further information send e-mail to: cerlpbr@tx.technion.ac.il. 
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Figure 1: Loci of fold and Hopf bifurcation in the parametric plane h, m. The three Hopf 
curves correspond to g = 0.4, 0.5 and 0.6. Solid grey lines show the supercritical, and 
dashed lines, the subcritical bifurcation. Outside the cusped region, the unique stationary 
state suffers oscillatory instability within the loop of the Hopf bifurcation locus. Within 
the cusped region, there are three stationary states, of which two, lying on the upper and 
lower folds of the solution manifold, may be stable. The solutions on the upper fold are 
unstable below the Hopf bifurcation line. 
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Figure 2: The stability boundary in the plane {go,yo), separating the regions of subcritical 
and supercritical Hopf bifurcation. The unphysical part of the curve below the boundary 
of positive frequency is shown by a dashed line. The solid circle at go = 0.5, yo — 2 marks 
the double zero point. 



